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ABSTRACT 

We analyze the physical processes of gap formation in an inviscid protoplanetary disk with an 
embedded protoplanet using two-dimensional local shearing-sheet model. Spiral density wave launched 
by the planet shocks and the angular momentum carried by the wave is transferred to the background 
flow. The exchange of the angular momentum can affect the mass flux in the vicinity of the planet 
to form an underdense region, or gap, around the planetary orbit. We first perform weakly non- 
linear analyses to show that the specific vorticity formed by shock dissipation of density wave can 
be a source of mass flux in the vicinity of the planet, and that the gap can be opened even for 
low-mass planets unless the migration of the planet is substantial. We then perform high resolution 
numerical simulations to check analytic consideration. By comparing the gap opening timescale and 
type I migration timescale, we propose a criterion for the formation of underdense region around the 
planetary orbit that is qualitatively different from previous studies. The minimum mass required for 
the planet to form a dip is twice as small as previous studies if we incorporate the standard values of 
type I migration timescale, but it can be much smaller if there is a location in the disk where type I 
migration is halted. 
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1. INTRODUCTION 



Disk-planet interaction is one of the important topics in the planet formation theory. A low mass planet embedded 
in a disk excites the density wave (Goldreich and Tremaine 1979), and the backreaction from the density wave causes 
the planet to migrate in the disk (e.g., Ward 1986, Tanaka et al. 2002). The excitation of density wave at Lindblad 
resonances can be understood by linear analyses, although it has recently been pointed out that non-linear effects are 
also important at corotation resonances (Paardekooper and Papaloizou 2009). For a high mass planet, the interaction 
between the planet and the disk becomes nonlinear, and the gap opens around the planetary orbit (Lin and Papaloizou 
1986ab, Ward and Hourigan 1989, Rafikov 2002b, Crida et al. 2006). 

Gap formation around the planet is important both theoretically and observationally. From the theoretical point 
of view, gap formation determines the regime of planetary migration. If there is no gap opening, planetary migration 
is in the regime called "type I" , where the interaction between the planet and the spiral density wave is important 
(Goldreich and Tremaine 1979, Ward 1986, 1997, Tanaka et al. 2002). Type I planetary migration timescale is 
considered to be faster than disk dispersal timescale, which poses a serious problem in the theory of planet formation. 
If the gap opens around the planet, the migration is in the regime called "type II", where the planet migrates as the 
disk accretion onto the central star occurs (Lin and Papaloizou 1986). The timescale of type II migration, which is of 
the order of viscous timescale, is generally longer than type I migration, and there may be a possibility for the planets 
to survive in the disk. From observational point of view, the gap around the planet may be able to be observed by direct 
imaging. Recent progress of disk observation by direct imaging has reached the stage that it is possible to compare the 
numerical simulations and observation directly (Mayama et al. 2009). Moreover, the dynamical interaction between 
the circumstellar dust or gas and the planet can be used to estimate the mass of a low-mass object embedded in the 
disk (Kalas et al. 2008). If the gap in the disk can extend to the disk scale, the gap structure can be a very good 
indicator of the existence of the planet. 

Conventionally, gap opening is understood as a balance between the torque exerted by the planet and the viscous 
torque (Lin and Papaloizou 1986b). In addition to the viscous torque, Crida et al. (2006) has shown that "pressure 
torque" also acts to balance the planetary torque. Gap opening processes are mainly investigated using one-dimensional 
model. Rafikov (2002b) calculate the evolution of disk surface density in the vicinity of the planet using a classical 
one-dimensional model (Lynden-Bell and Pringle 1974, Pringle 1981). He has shown that in the vicinity of the low- mass 
planet, the torque exerted by the planet is carried away by density wave0 , and the gap is not opened until the wave 
shocks to deposit the angular momentum to the mean flow. Considering the pressure effects as well as viscous effects, 
Crida et al. (2006) have derived the gap opening criterion generalized for arbitrary values of kinematic viscosity using 
two-dimensional numerical simulation (equation (15) of their paper). Using the criterion of Crida et al. (2006), the 
gap-opening mass for an inviscid disk reads 



where M p is the mass of the planet, M* is the mass of the central star, H is the scale height of the disk, and r is the 
orbital semi-major axis of the planet. Therefore, planets more massive than Saturn can open up the gap in the disk. 

However, recently, Li et al. (2009) have performed high resolution numerical simulations and showed that a partial 
gap is formed in the vicinity of the planet even if the mass of the plant is smaller than the mass limit given by equation 
([1} in case the disk viscosity is very low. This motivates us to investigate the physical processes of gap opening in an 
inviscid disk in detail. 

In this paper, we show that low mass planets can potentially open up a partial gap. We investigate the processes 
by means of analytical models taking into account weak non-linearity. We also perform numerical simulations to 
look at to what extent numerical calculations and analytic studies agree, and then finally we suggest the criterion 
for the gap opening in an inviscid disk. For analytic studies, we study the propagation of the spiral density wave 
and subsequent shock formation. We then investigate the mass flux in the vicinity of the planet by means of second- 
order perturbation theory. We note that the second-order perturbation is necessary to study the mass flux since it is 
essentially a second-order quantity. 

The plan of this paper is as follows. In Section [2j we describe the basic equations. In Section [3l we investigate the 
gap opening processes using a second-order perturbation theory. We show that the shock dissipation of density wave 
and the subsequent formation of specific vorticity can lead to the mass flux in the vicinity of the planet, resulting 
in the gap formation. We then show the results of numerical simulations in Section 2J In Section [5J we discuss the 
condition for gap opening in an inviscid disk. We also discuss that commonly used one-dimensional models of disk 
evolution may overlook gap opening processes in the vicinity of the planet. Section [B] is for summary. 



In this paper, we focus on the two-dimensional local shearing-sheet analysis for simplicity. We set up a local Cartesian 
coordinate system corotating with a planet. We take the origin of the coordinate system at the planet's location, and 
the x- and y-axes are the radial and the azimuthal direction, respectively. We use isothermal ideal hydrodynamic 
equations 




(1) 



2. BASIC EQUATIONS 




(2) 
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— +v-Vv = -— VS - 2Q p e z x v + 30^ - V^ p (3) 

where we have assumed that the gas is rotating at the Kepler velocity. Notations are as follows: £ is surface density, 
v is velocity field, c is sound speed, tt p is the angular velocity of the planet, and tp p is the gravitational potential of 
the planet. For ip p , we assume the form 

* - (4) 

where 67, M p , and e are the gravitational constant, the mass of the planet, and the softening parameter, respectively. 

Local shearing-sheet approximation is only an approximation of global model, and it may not be appropriate for 
investigating the global evolution of the disk structure. However, local shearing-sheet approximation and full global 
model share many essential physics in common. Excitation and the propagation of density wave can be understood 
using local approximation. We also show later that one-dimensional disk evolution model constructed from global 
model and local model are very similar. Local approximation also has the advantages that analyses are greatly 
simplified and high-resolution calculations become possible. 

Here, we comment on the two-dimensional approximation. It is known that three-dimensional processes are im- 
portant in disk-planet interaction especially when considering the immediate vicinity of the planet. For example, 
numerical simulations by Paardekooper and Mellema (2006) clearly show that disk structure around the planet is 
three-dimensional. However, for the study of density wave, which is the structure away from the location of the planet, 
the density perturbation by the planet is nearly two-dimensional. In this paper, we shall investigate the physical prop- 
erties of density wave excited by the planet in detail, and we discuss the gap formation processes that can be derived 
from the study of density wave. Therefore, we expect that essential physics can be captured by local two-dimensional 
model. 



3. ANALYTIC STUDY OF PLANETARY WAKE 

In this section, we investigate the disk-planet interaction and subsequent gap opening processes using analytic 
methods. We consider a weakly non-linear stage, where the mass of the planet is 

f^, (.) 

which is approximately smaller than the Saturn mass in case of Minimum Mass Solar Nebula (Hayashi et al. 1985) 
with H/r ~ 0.05 ■ We note that non-linearity of disk-planet interaction and gap opening is strongly related. From 
equation ([5]), the onset of the non-linearity is given by 

M >;>1.25x (6) 



M* ~ V 0.05 , 

where H = c/O p is the scale height of the disk. This non-linear criterion is analogous to equation ([1]), which is the 
gap opening criterion given by Crida et al. (2006) 
The overall picture of gap formation we suggest is as follows (see also Figure [T|) . 



1. Density wave excited by the planet will shock at some location away from the planet. 



2. The shock formation leads to the formation of specific vorticity. 



3. The change of specific vorticity results in net radial mass flux. This mass flux exists in the place closer to the 
planet, even at the place where the change of specific vorticity is not significant. 



4. Gap opens in the vicinity of the planet. 



We investigate the shock formation process and the mass flux separately. Shock formation processes are investigated by 
Goodman and Rafikov (2001) and Rafikov (2002a) using Burgers equation model, and the mass flux is investigated by 
Lubow (1990), using second-order perturbation theory. We show how these two theories can be combined to investigate 
the gap opening processes. 

Physically, it is natural that the gap opens when spiral density wave damps, since the angular momentum flux carried 
by the density wave should be transferred to the background flow. However, we shall show that the decay of the spiral 
density wave away from the planet can affect the mass flux in the immediate vicinity of the planet, in contrast to the 
model suggested by Rafikov (2002b), in which there is no gap formation in the vicinity of the planet in inviscid cases. 

3 Equation (01 is equivalent to the condition where Hill's radius becomes comparable to the disk scale height. Departure from linear 
calculations can be observed at this mass range, see e.g., Miyoshi et al. (1999) 
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3.1. Linear Analysis 

In this section, as a preparation for the non-linear analyses performed in the subsequent sections, we briefly summarize 
the results of linear study of the spiral density wave (e.g., Goldreich and Tremaine 1979, 1980). We denote background 
state by subscript "0", and the perturbation by 5: 



v 



T = T Q + 6T, 
3 



Sv 



Sv. 



(7) 

(8) 



The planet potential is regarded as a perturbation. The linearized equations are 

ST d_ 



dt 2 pX dy 



dx 
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(9) 
(10) 

(11) 



We now assume the stationary state in the frame corotating with the planet, d/dt 
y-direction. The perturbed values are given by 



5f(x,y)= ]T 5f(x)e 



ik y y 



0, and Fourier transform in the 

(12) 



n v £Z 



where Sf denotes ST, Sv x , or 5v y , and k y is the wave number in the y-direction. Assuming the periodicity in the 
y-direction, k y = 2irn y / L y , where n y is an integer and L v is the box size in the y-direction. The summation in the 
above equation is taken over the relative integers denoted by n y . It is possible to derive a single second-order ordinary 
differential equation (Artymowicz 1993). 



dx 2 



Sv v + 
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(13) 



and other perturbed quantities are given by 

ST _ 1 
% ~ D 

and 



Sv r = 

D 



where 
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~2~dx~ y 
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dx J 
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(14) 



(15) 



(16) 



The boundary condition for Equation (I13[) is that wave should propagate away from the planet. At the location 
away from the effective Lindblad resonance given by 



9 ^p k v_ x 2 



1 n l 

*2 - -£ = 0, 



4 c<- 

the approximate solution can be written analytically using WKB approximation . It takes the form, for |x| 

C{ky) 



5v v 



exp 



I . 3 Vipky 

±1 — 

4 c 



(17) 

00, 
(18) 



where upper sign is for x > and the lower sign is for x < 0, and C{k y ) denotes the amplitude (but not necessarily real) 
of the wave with mode k y . We note that since we consider the linear perturbation excited by the planet's gravitational 
potential at the origin of our coordinate system, the amplitude C(k y ) is proportional to the planet mass M p . From 
here on, using the symmetry of the shearing-sheet, we only consider x > without loss of generality. Using equations 
(fT4|) . (fT5|) , and (fT8|) , we can derive the asymptotic form for ST and Sv x at \x\ — > 00, 



ST 



D 



3n p 3 



exp 



3 i^p/Cf/ n 

l- x 

4 c 



(19) 



The exact solution of equation (1 1 3 I t is given by parabolic cylinder function, see Artymowicz (1993) for detail. 
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Sv x = c— , (20) 

where we have retained only the leading terms in x in equations (fT4")l and (fT5|) . Note that, in the real space, solutions 
are given by 

6f = Y / F(k y )x«exp l ~ ( 3x2X1 



(21) 



where H is the scale height H = c/fi p , g = 1/2 for <5£ and <5u x , g = —1/2 for Sv y , and F(k y ) is a function of k y . 
Therefore, along the line 

y = -~ + yo, (22) 

where yo is constant, perturbed quantities take the same values except for the dependence x q . Therefore, we can write 
the form of the WKB solution in real space, 

SE(x, y) = M p x^ 2 f(y + (3/4)x 2 /H), (23) 
Sv x (x,y) = M p x 1 / 2 g(y+(3/A)x 2 /H), (24) 

5v y {x, y) = Mpx-^hiy + (3/A)x 2 /H), (25) 

where we write the dependence on the planet mass explicitly, and /, g, and h are the functions that determines 
the form of the perturbation. We note that from equation (f20|) . / and g are proportional in the place where WKB 
approximation is valid. 

The exact linear solution can be obtained by solving equations ([9 |) -(|li p numerically with non-reflecting boundary 
conditions. The profiles of density and v x obtained in such a way are shown in Figure [5] In this figure, we calculate 
non-axisymmetric modes (k y ^ 0) and assumed GM p / He 2 = 1. The amplitude of perturbation is proportional to the 
planet mass. We note that the results obtained in equations ([23]) - ([25]) are based on the WKB approximation, and they 
are applicable only in the region \x/H\ > 1. 

3.2. Shock Formation 

Goodman and Rahkov (2001) performed a non-linear analysis of the propagation of the spiral density wave, and 
their local approach was later extended to the global model by Rahkov (2002a). Using several approximations based 
on the linear theory, they derived the Burgers' equation which describes the propagation of spiral density wave and 
concluded that the spiral density wave eventually shocks as it propagates in the radial direction. They derived that 

— 2/5 

the location of shock formation is proportional to M p . In this section, we derive this relationship using a slightly 
different consideration. 

Shock is formed when two characteristics cross. For two-dimensional supersonic steady flow, the gradient of charac- 
teristic curves is given by (Landau and Lifshitz 1959) 



dy\ _ v x v y ± c\Jv 2 - c 2 

dx ' ~ - 2 - 2 • [ b) 



± 



In the background state of the shearing-sheet, this gives 



<>!l \ . I ( -Kv 2 2 



For x > 0, (dy/dx) + is the perturbation propagating away from the origin. At x 3> (2/3)H, the characteristic curves 
are given by 

y~-~ + yo, (28) 

where j/o is a constant. It is to be noted that the outgoing characteristics coincides the curve of the same phase of the 
perturbation of the density wave (see equation (|22I0 . This is because the density wave is essentially the sound wave 
propagating on the disk. 

For the flow distorted by the perturbation of the planet, the characteristic curve is also distorted. Using the results 
of linear perturbation, the gradient of the characteristics is given by 



dy\ 1 



dx 



,.2 



Tc ( ^-nlx 2 - c 2 V + x J 5v x ± ''' )r " 



4 P ' J 2 P 1 ^0/^02^2 _ J\l/2 



((9/4)fi^2 



(29) 
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upto the lowest order of perturbation. Since the amplitude of Sv y decreases with x 1 I 2 , Sv y in the second term of the 
right hand side can be neglected for \x/H\ 3> 1. Assuming that the perturbed characteristic curve is given by the form 

3 x 2 

v{x) = -^jj+Vo + Sy{x), (30) 

where Sy(x) is given by 



3fl 

ax 

Approximating y(x) ~ — (3/2)x 2 /H + yo in the argument of 5v x , we have 



T- s y = ^-f xSv x(x,y(x)). (31) 
ax 2 c z 



5y{x) cx M p x 5 / 2 g(y ), (32) 

where g(y) is the function that appears in the linear solution of 5v x given by equation (1241) . From Figure [H g(y) has a 
zero point. In the vicinity of g(yo) = 0, g{y) is negative for y > yo and g{y) is positive for y < yo. This indicates that 
the separation between two characteristic curves shrinks compared to the unperturbed case. Assuming that when Sy 
exceeds a critical value, the shock forms, we obtain, for the location of the shock formation, 

x oc M~ 2/5 , (33) 

which is the same condition as derived by Goodman and Rafikov (2001). Since the flow is supersonic only at x > 
(2/3)H, condition 

2 



3* 



<x M~ 2 ' h (34) 



may be more appropriate. 



3.3. Second-Order Perturbation Analysis and Mass Flux 

In this section, we consider second-order linear perturbation in order to derive the mass flux in the vicinity of 
the planet. Although angular momentum flux can be calculated using the results of linear perturbation only, it is 
necessary to perform second-order analysis to derive the mass flux, since axisymmetric (k y = 0) mode of the second- 
order perturbation contributes to the mass flux. Lubow (1990) performed the time-dependent analysis using Fourier 
and Laplace transformation. In this paper, we calculate second order perturbation in real space. 

The equations for second-order perturbation are given by 

a_3 b\M» d_ (2) d_ (2) = 

dt 2 Up dy) So dx 6Vx + dy y ~ 

» f^^'f^,) (35) 



dx\ Y, x J dy \ S 



ov x gx dv x dv y dy dv x +c ^ ^ Eq (db) 



dt 2 ilp dy) 6 y +C dy S 



dVx dx dVy dVy d y dv y +c s dx s ' { ' 

The superscripts "(1)" and "(2)" denote the first- and second-order perturbation, respectively. We assume that the 
first-order results are already known. 
The mass flux is given by 

F M {t.x) = VoSvjP+SEMSvP, (38) 

where bar denotes the integral over y. Assuming the periodic boundary condition in the y-direction, it is possible to 
derive the equation for Tm which reads 

c 2 + ilpJ^M = S(t,x), (39) 



dt 2 dx 2 
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where S is the source term consisting of two parts, 



S(t,x) = S v (t,x) + —S t (t,x), 



where 



and 



S v (t, 



2fl r 



St(t,x) 



dt 



Eo 



(40) 
(41) 

(42) 



The term S v is related to the formation of specific vorticity. In two-dimensional ideal flow in a rotating frame, the 
specific vorticity conserves along the streamline 



d (V x v) z + 2Q T 
dt E 



= 0, 



(43) 



where d/dt = d t +v - V is the Lagrangian derivative. In the linear perturbation analyses, this reduces to d t — (3/2)il p d yt 
see also equations (|35|) -(|37 )) . The background value of the specific vorticity is f2 p /2Eo, and the linear perturbation is 



dt 2 p dy 
If there is no formation of specific vorticity, 



dy x dx y + 2 iZp E 



= 0. 



dy x dx y + 2 p So 



0. 



(44) 



(45) 



From equation (1411) . S v (t,x) becomes a total derivative with respect to y in this case, and therefore S v = 0. However, 
if there is a formation of vorticity by, for example, shock damping of the spiral density wave, this term can not be 
neglected. We also note that if stationary state is assumed a priori, and if there is no formation of specific vorticity, 
the source term S(t,x) is zero, leading to the zero mass flux (Lubow 1990, Muto and Inutsuka 2009a). 

However, if time-evolution effects are taken into account, we have non-zero mass flux. The solution for equation (|39|) 
is given by 

dt dx Q G(t,t Q ;x,x )S(t ,xo), (46) 



T M {t,x) 

c 

where G(t, x, xq) is the Green's function 

i J ^ ^c 2 (t - t ) 2 -(x- xo) 2 ) \x - x \ < c(t - t ) 
otherwise 



G(t,t ;x,x ) 



(47) 



where Jo is the Bessel function of zeroth order. From linear perturbation, it is possible to predict that the mass flux 
scales with M 2 , since the source term is the second order of perturbation. Later in this paper, we compare this result 
with numerical calculations to understand how much mass flux is excited by the planet. 
We now consider the model in which the source term is given by 







t < 



S{t,x) = 



So [Sd(x - x s ) - S D (x + x s )} t>0 



(48) 



where So > and x s > are positive constants, and Sd(x) denotes the Dirac's delta function. We later see that the 
source term is positive in the region x > 0, and negative for x < 0. This form of the source term is the simplest case 
where we can obtain an analytic solution for the mass flux. 
Changing the integration variable from (to,Xo) to (r, 8) via 



x — Xq 



^cHt-tof-i 



•'■ - ^o) 2 = — sin^, 



equation (1461) can be rewritten to 



F M {t,x) = 



dr I d8 sin 9 Ji 



S 



t , x — r cos ( 

c 



(49) 
(50) 

(51) 
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In order to see the mass flux only in the vicinity of the planet, we assume < x < x s . Substituting equation (|48[) . and 
integrate over r, we obtain 



' cos -1 [(x — x s ) fct] 
rcos~ 1 [(x-\-x s ) /ct] 

10 



cW 



tan 6* Jo ( X jj S tan 6* 



tan 9J U [ X ^J Cs t&nO 



(52) 



We take the limit t — > oo. Then, we can approximate cos [(x ± x s )/ct] ~ ir/2. Using the formula (Abramowitz and 
Stegun 1970) 

{•tt/2 



'0 JO 

where a > 0, the integration over 9 can be performed to obtain 



d9 tan 9 Jq (a tan 9) = I du- ~ Jo(au) — Kq(cl), (53) 

1 + u z 



2c z 



(54) 



where K is the modified Bessel function of the zeroth order. In the vicinity of the planet, x <C x s , the mass flux 
changes with Tm oc x, which indicates that the gap opens up. We note that divergence at x = ±x s is the artefact of 
our simplification where we have used delta function as a source term. 

3.4. Instability of a Disk with Surface Density Variation 

We have seen that if there is a source of specific vorticity, mass flux appears in the vicinity of the planet, and it 
leads to the change of surface density to open a gap. We now briefly discuss how this gap opening process ends. 

An inviscid disk with a rapid surface density variation is prone to a linear instability, which is referred to as Rossby 
wave instability (Li et al. 2000, de Val-Borro et al. 2007). In Appendix [Bl we show a brief outline of the linear stability 
analyses for a disk with a gap, and derive the necessary conditions for the instability. Gap induced by a planet in the 
disk naturally excites the variation of specific vorticity, and gap edges are likely places where such instability occurs. 
Rossby wave instability may stop the gap-opening processes described in previous subsections. We do not expect that 
the instability leads to the complete closing of the gap, since the planet always try to repel the fluid element by the 
excitation of the density wave. There may be at least "underdense region" around the planetary orbit, which may be 
called a "gap" . 

3.5. Summary of Analytic Study 

We have discussed the non-linear shock formation of density wave in Section 13.21 The formation of shock leads to 
the formation of specific vorticity in general, and the perturbation of specific vorticity can lead to the radial mass flux 
as discussed in Section 13.31 The gap formation may end by the onset of the linear instability of the disk with a gap. 
In order to investigate this qualitative picture of gap formation, we perform two-dimensional numerical calculations in 
the subsequent section. 

We note that in this picture of gap formation, the decay of angular momentum flux carried by the spiral density 
wave is not directly connected to the mass flux. For example, when we model the source of mass flux by <5-function 
in equation (j48|). we implicitly assume that the spiral density wave does not dissipate at |x| < x s , but we still end up 
with non-zero mass flux in this region, see equation (|54[) . One may consider that this result violates the conservation 
of angular momentum. However, it is actually satisfied, when time evolution effects are taken into account. We shall 
discuss further on the angular momentum conservation in Section 15.11 



4. NUMERICAL STUDY OF PLANETARY WAKE 
4.1. Numerical Setup 

We solve Euler equations ([2]) and ((3]) with isothermal equation of state using second-order Godunov scheme (e.g., 
Colella and Woodward 1984). In order to investigate the density wave propagation while resolving all the length 
scales (Hill radius, Bondi radius, and radial wavelength of the density wave), we need to use rather large box size with 
relatively high resolution. We choose the box size (L x ,L y ) = (16H,32H) with mesh number (N x ,N y ) = (512, 1024), 
which results in the resolution Ax = Ay = H/32. We have also used the box with different L y to investigate the box 
size effect. 

From equation (|21j) . the radial wavelength of the wave with mode k y may be approximated by 

The most important modes are k y H ~ O(l). For mode with k y H — 10 at x/H ~ 4, the wavelength of the density 
wave is given by A ~ H/4, and therefore, the radial wavelength is resolved by eight meshes if we use Ax = H/32. 
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Higher-order modes or the wave outside this regions are damped numerically. Therefore, we mainly use data within 
| at | < &H in the analyses of the numerical results. 

If we normalize the length by H = c/f2 p and the time by fi~ , the only dimensionless parameter in this calculation 
is the planet mass, or Bondi radius t-q/H = GM p /Hc 2 and the softening parameter e/H. We use the softening length 
used in previous local shearing-sheet calculations e = rjj/4, where r H is the Hill's radius 

?-H {M p \ 1/3 r p fGM p \ 1/3 (5g) 



H \3mJ H \3Hc 2 J 

(Miyoshi et al. 1999). We have performed calculations with five different planetary mass M p shown in Table [T] Note 
that Bondi radius, Hill radius, and the softening parameter are resolved (at least marginally) for the smallest mass 
model. We have increased the planet mass linearly from t£l p = to 12. The result does not depend significantly on 
this timescale if the timescale is longer than this. 

In the ^-direction, we use non-reflecting boundary used by FARGO (Barutcau 2008), modified for the shearing-sheet. 
For the y-direction, we adopt the periodic boundary condition. Shearing-sheet calculations performed by Miyoshi et 
al. (1999) or Tanigawa and Watanabe (2002) used a different boundary condition in the y-direction, which is the 
combination of Keplerian inflow and supersonic outflow. This boundary condition may be appropriate to study the 
flow structure only in the vicinity of the planet, but for the study of gap formation, this boundary condition is 
inappropriate because this forces the unperturbed gas flowing into the computational domain. We also note that 
periodic boundary condition in the y-direction is useful for the purpose of comparison with the linear analyses, which 
assume the periodicity in the y-direction. We have checked that our code reproduces the results of Miyoshi et al. 
(1999) well when the same parameter and the boundary conditions are used. 



4.2. Results 
4.2.1. Disk Structure and Evolution 

Figures [3] and 0] show the <5£ and 5v x at t£l p = 200 obtained by numerical simulations. For the numerical calculations 
with GM p /Hc 2 = 0.4, we can see that a gap is already formed. Figure [5] shows the evolution of the azimuthally- 
averaged density profile for GM p /Hc 2 — 0.1 and 0.4. It is possible to see that even for low-mass calculations, density 
gap is being gradually formed. In subsequent sections, we provide the physical interpretations of disk structure and 
evolution using analytical theories provided in Section [3] 



4.2.2. Spiral Shock Wave 

We first investigate the shock formation. Figure [6] shows the density profiles at various x at if2 p = 200 obtained by 
numerical calculations. It is possible to see the shock-like structure for the calculation with GM p /Hc 2 — 0.4, while 
for the run with GM p / He 2 = 0.1, the structure is not very obvious. 

It is possible to see whether dissipation acts or not by looking at the perturbation of specific vorticity. In the 
absence of dissipation, the specific vorticity conserves along the streamline, and since the background specific vorticity 
is constant in the calculation box, we expect that it is also constant in the presence of the planet. Note that any 
potential force, whether it is time-dependent or not, does not produce specific vorticity. Therefore, specific vorticity 
arises only if dissipative mechanisms come into play. In our numerical calculation, the dissipation is implemented in 
the shock-capturing scheme. If the formation of specific vorticity is driven by the formation of shock, relation given 
by equation ([34|) should be observed at least for the weak shock cases. 

Figure [7] shows the evolution of the perturbation of azimuthally averaged specific vorticity for calculations with 
GM p /Hc 2 = 0.1 and 0.4. It is possible to see that the specific vorticity is formed around x/H ~ 1.5. We later see 
that this specific vorticity becomes a dominant source of mass flux, which leads to the formation of a gap. 

For the numerical calculation with GM p /Hc 2 = 0.1, we also see the formation of specific vorticity around x ~ 0. 
This is because the fluid elements in the vicinity of the planet orbit around it, and as a result, a vortex is formed 
just around the planet. The formation of a vortex at the planet location causes the vortex circulating in the opposite 
direction, which makes a horseshoe orbit in the corotation region. For the calculation with GM p /Hc 2 — 0.4, similar 
thing happens, but the vortices formed by shock dissipation are stronger. In subsequent sections, we show that vortices 
just in the vicinity of the planet location does not produce a significant amount of mass flux. 

Figure [8] shows the azimuthally- integrated specific vorticity perturbation as a function of {x- (2/3)H)M p /5 . If shock 
dissipation occurs as equation (|34[) . the peak of the specific vorticity perturbation comes at the same location of the 
horizontal axis. It is possible to observe that the for calculations with GM p /Hc 2 < 0.2, equation ([34|) is marginally 
satisfied. For calculations with larger planet mass, it is not the case. This is because the shock formation occurs 
immediately after the excitation of the wave at x ~ (2/3)H, and therefore in the unit of (x — (2/3)H)M p ^ 5 , the shock 
occurs relatively further away. 



4.2.3. Mass Flux 

We now look at the radial mass flux T,v x excited by the planet. In figure [9j we show the azimuthally averaged mass 
flux obtained by numerical calculation. It is possible to see that there is a spike of the mass flux just in the vicinity of 
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the planet, and then the mass flux depends linearly with the distance from the planet within \x/H \ < 2. The spike in 
the vicinity of the planet is due to the gas falling onto the planet. This causes the spike of the surface density profile 
in the vicinity of the planet as shown in figure [5] We expect that this feature depends on the treatment of the planet, 
and such mass flux should be investigated carefully when one wants to look at the processes such as gas accretion onto 
the planet. However, as we will see later, the feature in this spike region does not affect the mass flux outside this 
region. In the discussion below, we focus on the structure outside this spike region. 

We first compare the results obtained by numerical simulation and analytic calculation. Figure [10] compares the 
mass flux derived by numerical calculation and by equation (|46[) . Since it is not possible to predict the amount of 
specific vorticity only from linear perturbation theory, we have used the results of numerical calculation in obtaining 
the source term S(t,x). It is possible to observe that the second-order perturbation theory is valid for calculations 
with low-mass planet, while for calculations with GM p j He 2 — 0.4, the perturbation theory fails to explain the amount 
of mass flux. We find that the second-order perturbation theory can be used for GM p /Hc 2 <, 0.2. 

We further investigate how this mass flux is formed. We first look at the source term of the mass flux given by 
equations (|4Tj) and (|42j) . We have used the results of numerical calculations to calculate the perturbed values that 
appear in these two equations. Figure [TT] shows the evolution of the source of mass flux given by fill) and (|42l) for 
the calculation with GM p /Hc 2 = 0.1. We see that the term with dfS(t,x) is not significant after the planet mass is 
fully introduced at tfl p = 12. The term S v (t,x) depends on time strongly in the vicinity of the planet, while nearly 
time-independent contribution is observed from the region \x/ H\ > 1.5. The source arising in the vicinity of the planet 
comes from the vortex making a horseshoe orbit, while the source at \x/ H\ > 1.5 comes from shock formation. 

In order to look at which contribution excites the mass flux more effectively, we calculate the mass flux by artificially 
cutting the source term at \x/H\ > 1. Figure [T"2l compares the mass flux obtained by doing so and the mass flux 
obtained by using the full source term. We have used equation PS)) with the source terms obtained by numerical 
calculation to derive the mass flux. It is clear that the source in the vicinity of the planet does not contribute to the 
mass flux very much. Therefore, we conclude that the mass flux arises from the specific vorticity generated by the 
shock dissipation of density wave. 

We now discuss the dependence of mass flux on the box size in the y-direction, L y . Figure IT3l shows the azimuthally 
integrated mass flux J dyT,v x with GM p /Hc 2 = 0.2 for different box length L y . These mass fluxes are directly 
calculated from the numerical simulations. The lines show a perfect match for different box sizes. This can be 
explained if we notice that the angular momentum flux carried by the density wave does not depend on the box size. 
As the spiral density wave damps, the angular momentum carried by the wake is deposited to the background disk to 
drive the mass flux. The total amount of the angular momentum deposited over the y-direction does not depend on 
the box size, since the amount of angular momentum flux carried by the wave does not depend on the box size. 

We now discuss the timescale for the gap-opening in the vicinity of the planet. As stated at the end of Section l3~3l 
the mass flux in the vicinity of the planet is expected to be proportional to x, and if second-order analysis is valid, 
we expect that the mass flux is proportional to the square of the planet mass. We therefore fit the mass flux within 
\x/H\ < 2 using the function of the form 

where K is a constant that is derived by fitting the results of the numerical simulations. We found K = 4.4 x 10~ 3 
can fit the numerical results upto GM p /Hc 2 <J 0.4 within 15% error for \x/H\ < 2. Figure [§] compares the mass flux 
obtained by numerical simulation and the fitting function (|57[) , 
Since the gap opening timescale may be determined by 

™ (l/L tf )/d tf E ' l " 8j 

and the denominator is approximated by So; the gap-opening timescale is proportional to L y . Using the result of 
equation (1571) . we have 

Note that this timescale applies to the initial phase of gap opening, and it does not state about how deep the gap will 
be. The box size in the y-direction may be interpreted as a circumference of the disk. The above expression can be 
rewritten as 



4.2.4. Instability of Gap Profile 

In Section 13.41 we have seen that the disk with a density bump can be prone to linear instability. Figure [14] shows 
the evolution of azimuthally averaged density perturbation and mass flux in the calculation with GM p /Hc 2 = 0.6. 
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We see that the mass flux changes the sign very rapidly, and the gap depth saturates when approximately half of the 
original mass is depleted. The width of the gap we have obtained is of the order of several scale height of the disk. 

It seems that the gap width still increases gradually. This may be explained quantitatively as follows. The perturba- 
tion by the planet tries to deposit the angular momentum to the disk in such a way that the fluid elements to migrate 
away from the planet. Turbulence may try to redistribute the angular momentum, but at the same time, diffusion 
in the radial direction occurs because of the turbulence. Therefore, not all the angular momentum deposited to the 
disk is absorbed by turbulence, and a part of the angular momentum goes to the background, and the gap gradually 
becomes wider and wider. It may be possible to form a wider gap. The final state of the gap may not be able to 
capture in the local shearing-sheet calculation since the radial extent of the calculation box is not very wide. 

Li et al. (2009) obtained a partial gap (approximately half of the gas depleted) with the width of about 5 — 7H for 
the parameter GM p /Hc 2 ~ 0.7 (see Figure 2 of their paper). Although it may be possible that the final width and 
the shape of the gap is not well-determined in the shearing-sheet calculations, the gap formation process described in 
the previous sections may be valid. We interpret the gap depth and width obtained in our calculation as minimum 
values. Planets may potentially be able to open a deeper and wider gap. We note that Crida et al. (2006) and de 
Val-Borro et al. (2007) obtained deeper and wider gap for a Jupiter mass planet, although they assumed very low or 
zero viscosity. 

5. DISCUSSION: GAP FORMATION IN A PROTOPLANETARY DISK 

In this section, using the results of numerical calculations, we discuss the validity of conventional one-dimensional 
model for gap formation and the condition of gap opening. 

5.1. One- Dimensional Model 

The gap formation processes are modeled using one-dimensional disk evolution model. Thommes et al. (2008) 
uses a model in which the torque exerted at the Lindblad resonances are directly deposited to the disk to change 
the semi-major axis of the fluid particle. Rafikov (2002b) uses a one-dimensional model based on Lynden-Bell and 
Pringle (1974) type surface density evolution equation to model gap formation processes. We investigate whether such 
treatment can reproduce the evolution of the density profile obtained in two-dimensional calculations. 

Although we have performed shearing-sheet analyses, we first note that global calculation and local shearing-sheet 
calculation share very similar property if we consider the azimuthally-averaged one-dimensional model. This can be 
investigated by comparing the one-dimensional model derived usi ng glo bal m odel and local model (see Appendix [A| . 

Exact one-dimensio nal e volution model is given by equations (|A12|) and (|A15j) in local approximation. Whether 
the evolution model (|A17[) can be used depends on whether the approximation made in deriving this equation is 
appropriate. This can be seen whether mass flux, angular momentum, and the torque exerted by the planet are 
related as equation (|A16|) . In Figure [15] we compare each term of equation (|A16[) . We observe that the relationship 
of this equation is not sa tisfied, indicating the time evolution of T,5v y occurs. We have checked that the exact time 
evolution equation (|A15|) is satisfied. 

We therefore argue that time evolution of density and rotation profile is given by equation (|A12|) and (|A15|) , and in 
constructing one-dimensional model, it is necessary to specify the model of the mass flux, as well as the torque and 
angular momentum flux. It is especially important in considering the gap formation in the vicinity of the planet orbit. 
For example, Rafikov (2002b) predicts nothing happens in the vicinity of the planet in an inviscid model unless the 
spiral density wave shocks. In two-dimensional simulation, in contrast, we also observe gap opening even in low-mass 
planet calculations. As the gap opens, rotation velocity of the gas is also modified in order to balance the pressure 
gradient, as Crida et al. (2006) pointed out. 

5.2. Gap Opening Condition in an Inviscid Disk 

We have seen that low-mass planets can potentially open a gap, or at least a ring with low surface density, in an 
inviscid disk. We now discuss the minimum mass of the planet that can open the gap. Low mass planets are prone 
to type I migration and therefore, if the planet migrates before it opens a gap, gap formation is impossible. Since 
the width of the gap is the order of scale height H , we compare the timescale of the planet migration over length H 
and the gap formation timescale. The planet migration timescale in an isothermal disk is estimated by Tanaka et al. 
(2001) as 

_i _„ M n 5>2 /rnX 3 



rZt ~ 3ft 



mig 



p^^PO . (61) 



Note that the power of H/r p is 3 because we consider the migration over the radial length H . For gap formation 
timescale, we use the results of numerical calculations given in Section 14.21 The gap opening timescale T~Jp is given 
by equation (|59")l . If T m ^ g < T~^ pl we expect the gap will open in the vicinity of the planet. Comparing T m j g and r gap , 
we obtain 

^ > 2.o x io- f^V (hdn\ J (J^) 2 (MA 1 . (62) 

M* V 0.05 J \ 32 J 2 x 10 3 gcm- 3 UAU7 \M Q J V 1 

We have used typical values of protoplanetary disk at 1AU to estimate the number. We note that the gap-opening mass 
scales with the box size L y . If we assume L y = 2irr p and H/r p = 0.05, L y /H = 125.66, which gives M p /M* ^ 8x 10~ 5 , 
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which is below the criterion given in previous studies, e.g., equation ([T}. Equation (|62[) rewritten by using 27rr p reads 
(see also eauation (|60)) ') 



It is interesting to compare equation (|62p with the previous criterion (|TJ) . Equation (I62[) is derived by comparing 
the gap opening timescale and type I migration timescale, and this is qualitatively different from the way the previous 
criterion is derived. For example, the gap opening criterion derived in this study, equation (|62[) . depends on the disk 
mass, since type I migration timescale scales with the disk mass. We also note that the type I migration timescale 
may be much longer if the corotation torque and non-barotropic effects are considered (Paardekooper et al. 2010ab). 
In this case, our results indicate that the gap may be opened for lower mass planets, although more thorough analyses 
of non-barotropic effects on spiral density wave and mass flux are necessary. 

If the disk is in turbulent state, the effective viscosity can also act to fill the gap. If we use standard a-prescription 
for turbulent viscosity v = acH, the timescale of viscous diffusion over length scale ~ H is 

r-J ~ att p . (64) 

If T v is is shorter than r ga p, we expect that the gap is filled. This leads to another condition, 

a* 1.1 x 10- f^V 6 ( hdEY 1 , (65 ) 

V 0-05 J V 2 x !0~ 5 / V 32 J ' v ' 

or if we use 27rr p instead of L y , 

a < 2.8 x 10- (Wf7 .^Vf ( 2 -p). (66) 

V 0.05 / V 2 x !0~ 5 J \ L y J v ' 

The disk must be quiet in order for a low-mass planets to open up the gap, but we note that the condition is very 
sensitive to the disk aspect ratio. We note that the previous results for gap-opening criterion derived by Crida et 
al. (2006) is somewhat different from equation (|65p. Since we derive equation (|65[) using a crude order-of- magnitude 
estimate for viscous diffusion, it is necessary to perform more systematic parameter study in order to investigate how 
much viscosity is necessary to halt gap-opening by a planet. 

The gap-opening criteria, equations (|62[) and (|65p . are more exactly the conditions for the "formation of under-dense 
region around the planetary orbit" . We compare the timescale for the decrease of surface density, which is derived 
from the mass flux induced by the planet, and the type I migration timescale or viscous diffusion timescale. Therefore, 
these conditions concern the initial phase of the gap-opening, and the timescale argument does not concern how deep 
the gap will be as a result of long-term evolution. The condition by Crida et al. (2006) is, on the other hand, the 
condition for the formation of "deep gap" , or more quantitatively the formation of the gap where 90% of the original 
amount of gas is depleted. Therefore, the conditions are not necessarily in contradiction with each other. For example, 
in the case of the numerical calculation with GM p /Hc 2 = 0.6, 3H/4rn ~ 1.3, which is slightly above the gap-opening 
criterion by Crida et al. (2006). In figure [14j we observe a gap where approximately 50 — 60% of the original gas 
mass is depleted, which is just below the depth of the gap which is predicted by Crida et al. (2006) This may suggest 
that the prediction by Crida et al. (2006) and our calculation is in qualitative agreement, although the setting of our 
calculation is different from theirs. Our condition, however, still has a new aspect since it states that the gap-opening 
process may be different in case of inviscid case. Non-linear evolution of density wave and the feedback onto the disk 
can be important when we consider an inviscid disk. 

6. SUMMARY AND FUTURE PROSPECTS 

In this paper, we have investigated the disk-planet interaction and the evolution of surface density profile of the disk 
using analytic methods and high-resolution numerical calculations within the framework of shearing-sheet approxima- 
tion. We have shown some analytic framework to understand the non-linearity of disk-planet interaction. Formation 
of specific vorticity by shock dissipation of density wave can be a source of disk mass flux in the vicinity of the planet, 
which results in the formation of gap around the embedded planet. We estimate the conditions of the formation of the 
underdense region, or gap formation, around the planetary orbit, and it is indicated that a planet with 20-30 Earth 
mass at 1AU will open a gap in a standard Minimum Mass Solar Nebula with H/r ~ 0.05. We note that our condition 
applies to the initial phase of the gap opening, and the condition depends on the type I migration timescale, which 
can be much longer than that given by Tanaka et al. (2002) formula. If type I migration is halted, much lower mass 
planets can open a gap. In an inviscid case, which we have explored in this paper, the width of the gap is of the order 
of the disk scale height. We note that we have not discussed the final gap depth in detail in this paper, although it is 
indicated that the gap opening process ends by the onset of hydrodynamic instability when approximately half of the 
mass is depleted. One may call such a shallow underdense region a "dip" or "partial gap" , rather than gap. We have 
also discussed that classical one-dimensional model of disk evolution fails to describe gap opening process, and more 
refined one-dimensional is necessary to explain gap opening. 

We have focused on the two-dimensional local shearing-sheet calculations. We have discussed that the final width of 
the gap may be larger than we have obtained in this calculation. The gap depth and width we have obtained may be 
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interpreted as minimum values. Global calculations with high resolution are necessary to construct more quantitative 
model for gap formation. We also note that the discussion using the specific vorticity we have mentioned a number of 
times in this paper may be appropriate only in the two-dimensional calculations. Although we expect that the two- 
dimensional model is adequate for the investigation of spiral density wave qualitatively, three-dimensional calculations 
might result in the quantitatively different condition for gap opening. 
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APPENDIX 

DERIVATION OF ONE-DIMENSIONAL MODEL 

In this section, we compare one-dimensional dynamical evolution model derived from global and local model, and 
show that they share very similar properties. 



Global Model 

We use a cylindrical coordinate system (r, <fi) . In an inertial frame, a set of fluid equations in two-dimensional global 
model is given by 

f + ll^) + l^)=0, (Al) 

dv r | v dv r | v$ dv r v l _ I dp di/j p 
dt dr r d<j) r S dr dr dr 

dv^ ^ ^ dv^ | dv r v r v^ I dp 1 di/i p 

dt dr r dip r Er d(j) r d(f> ' 

where is the gravitational potential of the central star which is assumed to be axisymmetric. The rotation profile 
of the background disk is given by 

r n 2 ( r ) = (A4) 
dr 

We now derive a one-dimensional model for disk evolution by taking the azimuthal average of equation of continuity 
(|Alj) and conservation of angular momentum, which is essentially equation (|A3j) . We decompose the azimuthal velocity 
into background part and perturbation, 

v$ = rVL(r) + Svj, (A5) 

but we do not use linear approximation. 
Azimuthally averaged (denoted by bar) equation of continuity is 

§ + i^(r^)=0, (A6) 
at r or 



and the azimuthally averaged angular momentum conservation is 
where 



3 (r£ U0 ) + ~ (r 3 ftSX, + r^VrSv^) = T(r,t), (A7) 



T(r,t) = -£d^ P (A8) 

is the torque exerted by the planet. 
Decomposing as equation (IA5|) . equation (|A7[) can be further rewritten, with the aid of equation (|A6|) . 



^ (rm^) + ~ (r^VrSv^) = ~(r 2 n)'T(r, t), (A9) 

where / denotes the derivative with respect to r. 

Upto equation (|A9[) . our treatment is exact. In the standard one-dimensional model, approximation rfl ^> Sv r is 
used to neglect the time-derivative of equation (|A9[) to obtain the relationship between the mass flux and angular 
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momentum flux (Balbus and Papaloizou 1999). This approximation means that one neglects the evolution of the 
rotation profile. We then obtain the relation between mass flux, angular momentum flux, and torque, 



ld_ 

r dr 



(EtV^tv,) + T(r,t), 



and the evolution of surface density is obtained from equation (|A6 



9E 
dt 



ld_ 

r dr 



1 







(r 2 n)' dr 



1 d 


~rT(r,t)~ 


r dr 


(r 2 ny _ 



(A10) 



(All) 



If a-prescription is used for the angular momentum flux and the planet is absent, this is the model by Lynden-Bell 
and Pringle (1974). Rafikov (2002b) used this equation to study gap formation. 

Local Model 

It is possible to derive equations analogous to global model in averaging over the y-direction in the shearing- sheet 
approximation. Taking the y-average of equation of continuity, we obtain 



dZ d — 







and from the equation of motion in the y-direction in a conservation form, 

d 



dt^ 



3 p dx 



— T,V X 5Vy 



Ti oc (t,x) 



where 



T loc = -T,d y ip P - 

Decomposing the first term in equation (|A13[) into the background and perturbation, we obtain 

di + W 



— Y.5v y + —Hv x 8v y = --Q p Y,v x + T loc (x). 
This corresponds to equation (|A9I) in the global model. 



If we approximate that |(3/2)f2 p x| ^> 5v y , mass flux and angular momentum flux are related by 

dx' 

we obtain 



-VL v Y,v x = -—T,v x 5v y + T loc (x). 



as d_ 

dt dx 



^d-x^ Vx5Vy 



(A12) 

(A13) 
(A14) 

(A15) 

(A16) 
(A17) 



which corresponds to equation (|A11[) in global model. 
We also note that the equation for mechanical energy loss is also analogous in global and local model. 



LINEAR STABILITY ANALYSIS OF A DISK WITH A GAP 



In this section, we show the outline of the linear stability analysis of a disk when the surface density structure is not 
uniform. 

We consider a disk without a planet for simplicity. The equations we consider are equations @ and without 
We assume that the background disk is axisymmetric with density profile Eo(x). We denote the background values 
with subscript "0" . In the background state, the pressure gradient must be balanced by Coriolis force and therefore, 



and 



Vyi0{x) = --n pX+ — jr^^UW, (Bl) 

v x , Q - 0. (B2) 

We now consider linear perturbation. Perturbed values are denoted by S and consider the solution proportional to 
exp[— icut + ik y y}. Linear perturbation is then 



% + %^x- 5Vx + d^ SVx+tk ^ = °' 



2 d <5£ 

iujov x + c — — — 2\l D dv v 

dxT, 1 y 



0. 



(B3) 
(B4) 
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where uj(x) = w — k y U{x). From these equations, we can derive a single second-order ordinary differential equation 
for $ ee Eo<5w x , 



d_ 

dx 



1 



Sq (a) 2 — c 2 fc 2 ) dx 



feu 



where 



20 p £ \Y> {Cj 2 - c 2 k 2 ) ) c 2 S (£ 2 - c 2 A; 2 ) 



k 2 (x) = 2f2p(2r2 p + C/'(.t)) 



$ = 0, 



(B6) 
(B7) 



and / denotes the derivative with respect to x. Equation (|B6() with an appropriate boundary condition constructs an 
eigenvalue problem for oj. Since full analyses of equation (|B6I) is not the scope of this paper, we briefly outline the 
qualitative results about the instability of gap pr ofile that can be derived from equation (|B6|) . 
For an axisymmetric mode (k y = 0), equation (|B6|) becomes 



d_ 

dx 



1 d$ 

So dx 



2 2 

Ul — K 



= 0. 



(B8) 



We multiply $* to this equation and integrate over x. Assuming that the perturbation vanishes at \x\ — > oo and 
integrating by part, we obtain 



/ 



dx^- 



d$ 



dx 



dx- 



i*r = / dx 



1*1 



(B9) 



We first note that from this equation, the eigenvalue uj 2 must be real. For instability, w 2 must be negative and 
therefore, k 2 must be negative at some point in x. This necessary condition for instability is the well-known Rayleigh 
criterion. 

We can further proceed by using the independent variable 



1 



-A>. 



(B10) 



Equation (|B8[) now becomes 



Defining 



and 



dx 2 



+ 



V(x) 



n 2 
p 



i 



c 2 V^dx 2 v l " 

equation (|B11[) can be looked as a Scrodinger equation with energy E and potential V(x), 

^ + [E-V(x)}g(x) = 0. 



(Bll) 

(B12) 
(B13) 

(B14) 



If there is a "bound state" with energy eigenvalue E < 0, the system is unstable. Necessary condition for instability 
is therefore that there exists a point where V(x) < 0, or 
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,1/2 dx 2 



sy 2 <-^ 



n 2 
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(B15) 



Therefore, if there exists a sharp pressure bump, the system is unstable against axisymmetric perturbation. 

For a non-axisymmetric mode (k y ^ 0), the analysis becomes more complex. However, we can make a progress when 
we assume that there is a mode confined close to corotation, Co = 0. Assuming uj ~ so u> — c 2 k 2 ~ —c 2 k 2 , equation 
(1B6|) becomes 
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Multiplying $* to this equation and integrate over x, we obtain 



dx 



1 



dx 



c 2 kl + k 2 



Taking the imaginary part of this equation, we have 

lm(uj) / dx 



1 ki, 
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(B17) 



(B18) 



Therefore, should imaginary part of w exist, (k 2 /So)' must change the sign somewhere. This is Rossby wave instability 
previously investigated by a number of authors (Lovelace and Hohlfeld 1978, Lovelace et al. 1999, Li et al. 2000). Our 
analysis presented here for non-axisymmetric disk is the shearing-sheet version of that given by Lovelace and Hohfeldt 
(1978), although we have used a different variables from their analysis. We note that k 2 /T, is proportional to the 
background vortensity, which is (2f2 p + t/')/So in the shearing-sheet approximation. 
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1 . density wave shocks 




Fig. 1. — Overall picture of the physical processes involved in gap opening. The density wave launched by the planet eventually shocks 
as it propagates in the radial direction. Dissipation of the density wave then leads to the formation of the specific vorticity. The specific 
vorticity leads to the mass flux around the planet, which results in the gap formation. 
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Fig. 2. — Linear density perturbation i5E/So (left) and Sv x (right). We use GM P /Hc 2 = 1, but the perturbation scales with the mass of 
the planet. 
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Fig. 8. — Perturbation of specific vorticity (in arbitrary unit) as a function of [x — (2H/3)) X M p 
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result, GM /He =0.1 
fitting function ■ 



lumerical result, GM /He =0.4 +■ 
fitting function 



FlG. 9. — Comparison between the mass flux obtained from numerical calculation (dots) and the fitting function J57} with K = 4.4 x 10 3 . 
We plot the value (GM p f Hc 2 )~ 2 (1/L y ) f dyTiV x . Left panel shows the results with GM P /Hc 2 = 0.1 and the right panel shows the results 
with GM P /Hc 2 = 0.4. The horizontal axis shows x/H. 
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Fig. 10. — Comparison of the mass flux obtained by numerical calculation and second-order perturbation theory. Left panel shows the 
result with GM P /Hc 2 = 0.1, and the right panel shows the results with GM P /Hc 2 = 0.4. 



Gap Opening in Inviscid Disk 



27 




Fig. 11. — The evol ution of the source terms of mass flux obtained by the calculations wit h GM P /Hc 2 = 0.1. Left panel shows S v (t,x) 
given by equation 14H . Right panel shows dtSt(t,x) where St(t,x) is given by equation 14_'l . 
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Fig. 12. — Mass flux obtained by second order perturbation theory. Solid line shows the case where the source term is restricted to 
\x/H\ < 1, while the dashed line shows the result with full source term. 




Fig. 13. — Integrated mass flux J dyT,v x for box sizes with L y /H = 16,32,64, 128. The mass of the planet is GM p / He 2 = 0.2 and the 
data at tQ p = 100 are used. 
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-0.001 1 ' ' ' ' ' 1 

-4 -2 2 4 

x/H 

Fig. 15. — Comparison of e ach te rm in equation IIA16I ). Data with GM p / He 2 = 0.2 at ttt p = 100 is used. The line "(mass fiux)/2" is 
the left hand side of equation (I A 161 ) . "derivative of angular momentum flux" is the first term of the right hand side (sign inverted), and 
"torque" is the second ter m of t he right hand side. The line "total" is the value when all the terms are taken to the left hand side, and 
should be zero if equation ||A16|) is strictly satisfied. 
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TABLE 1 

Mass parameters used in numerical calculations 



Model Number 


GMp/Hc z 


Planet Mass a 


r H /H 


1 


0.05 


1.875M® 


0.26 


2 


0.1 


3.75M e 


0.32 


3 


0.2 


7.5M© 


0.41 


4 


0.4 


15M e 


0.51 


5 


0.6 


22.5M ffi 


0.58 



a Assuming 1AU and H/r p = 0.05. 



